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We study two continuous variable systems (or two harmonic oscillators) and investigate their 
entanglement evolution under the influence of non-Markovian thermal environments. The continuous 
variable systems could be two modes of electromagnetic helds or two nanomechanical oscillators in 
the quantum domain. We use quantum open system method to derive the non-Markovian master 
equations of the reduced density matrix for two different but related models of the continuous 
variable systems. The two models both consist of two interacting harmonic oscillators. In model 
A, each of the two oscillators is coupled to its own independent thermal reservoir, while in model 
B the two oscillators are coupled to a common reservoir. To quantify the degrees of entanglement 
for the bipartite continuous variable systems in Gaussian states, logarithmic negativity is used. We 
find that the dynamics of the quantum entanglement is sensitive to the initial states, the oscillator- 
oscillator interaction, the oscillator-environment interaction and the coupling to a common bath or 
to different, independent baths. 
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I. INTRODUCTION 

Quantum information and computation is a nascent 
and interdisciplinary field that exploits quantum effects 
to compute and process information in ways that are 
faster or more efficient than or even impossible on con- 
ventional computers or information processing devices. 
This field has typically concerned itself with the manipu- 
lations of discrete systems such as quantum bits (qubits). 
Recently, the extension to continue variables such as po- 
sition, momentum, or the quadrature amplitudes of elec- 
tromagnetic fields has led to the illuminating concept of 
continuous variable quantum information processing [lj. 
This includes the experimental realization of quantum 
teleportation 0, [IJland the demonstration of quantum 
key distribution p, [g] for continuous optical fields, and 
the successful definition of the notion of universal quan- 
tum computation over continuous variables 

Advances in current technology have allowed the fab- 
rication of very small mechanical cantilevers or oscil- 
lators with high frequencies, and allowed their opera- 
tion and manipulation at very low temperatures 0, Q. 
In the regime when the individual mechanical vibra- 
tion quanta are of the order of the thermal energy, 
the motion of the mechanical oscillators is close to or 
on the verge of the quantum limit. Recently inves- 
tigations devoted to observing quantum effect in the 
truly solid-state mechanical oscillators have been re- 
ported 0, M HH, Ei EI El EE M Ei El M- Th e 

Hilbert space of quantized electromagnetic fields is equiv- 
alent to the Hilbert space of the quantum Harmonic os- 
cillators. Thus, in addition to quantum optics system, it 
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may be possible to implement continuous variable quan- 
tum information processing in the nanomechanical oscil- 
lator systems. This is particularly interesting in that it 
provides a stepping stone towards quantum state control 
and a platform to explore the transition from quantum 
to classical world in mechanical systems that consist of 
many-million atoms [H El, El, El El, El • 

Quantum entanglement has been considered as key re- 
sources in many aspects and applications of quantum in- 
formation processing. In the real world, quantum co- 
herence and entanglement of quantum systems will in- 
evitably be influenced and degraded by external environ- 
ments. There have been several investigations of decoher- 
ence and quantum entanglement of continuous variable 
systems under open s y stem d y namics in the literature 

©llilS|2l&|2a|2S|2llMl2ll33,p. But in 

those investigations, Markovian approximation or (and) 
rotating-wave approximation (RWA) is (are) assumed. 
However, if the short time interval or regime, comparable 
with the environment correlation time, is concerned, or 
if the environments is structured with a particular spec- 
tral density, then the non-Markovian environmental ef- 
fect could become significant. For example, in the case 
when a high-speed quantum information processing is 
required, the non-Markovian effect becomes important 
since the typical characteristic time of the relevant sys- 
tem may be comparable with the reservoir correlation 
time. Besides, when the typical system characteristic 
time is comparable with the decoherence and dissipation 
times, the other approximation, RWA, widely used in 
quantum optics master equation [32| to describe open 
quantum system may not apply. This is particularly the 
case for nanomechanical oscillators (or beams) as their 
fundamental vibration frequency Q could currently just 
reach a few GHz [33[ , still much smaller than the optical 
frequency of 10 15 Hz. Thus a detailed investigation of 
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FIG. 1: (Color online). The schematic illustration of the two 
models investigated. The two models both consist of two 
interacting harmonic oscillators. In model A, each of the two 
oscillators is coupled to its own independent thermal reservoir, 
while in model B the two oscillators are coupled to a common 
reservoir. 

the non-Markovian quantum entanglement dynamics in 
a more general setting without RWA is demanding. The 
main purpose of this paper is to present such a detailed 
analysis. The analysis is, as mentioned, of great impor- 
tance to quantum nanomechanical oscillator systems and 
its relevance to quantum optics systems is also obvious. 

We study, in this paper, two harmonic oscillators in 
the quantum domain and investigate their entanglement 
evolution under the influence of thermal environments 
[HI GJ, [M HE H3, HI. Two different but related mod- 
els of harmonic oscillator systems are investigated (see 
Fig. [TJ. Model A consists of two interacting harmonic 
oscillators, each coupled to its own independent reser- 
voir. The two oscillators may be envisaged to be suffi- 
ciently spatially apart and may thus be relevant to the 
related setup for applications in quantum communication 
and teleportation. Model B also consists of two interact- 
ing harmonic oscillators, but both coupled to a common 
reservoir. They may be spatially close and could be use- 
ful for possible applications in quantum computation or 
other quantum information processing tasks. The cou- 
pling between the two oscillators, and the coupling be- 
tween the oscillators and the environments are, if present, 
all bilinear in their respective positions (or coordinates). 
In practice, there may not be direct interaction between 
two optical fields, and this can be simply achieved by set- 
ting the coupling between the two oscillators to be zero 
in our models. On the other hand, a controllable and 
tunable interaction between two nanomechanical oscilla- 
tors can be introduced by applying a voltage made from 
the metallic film fabricated on their surfaces [lfj, [l9|, [3J] . 
Thus the two models discussed here are applicable to 
quantum nanomechanical oscillator systems and quan- 
tum optics systems. 

In the context of two modes of electromagnetic field 
embedded in a thermal environment, Ref. [25j derived a 
condition which states that if the state of the two modes 
is initially sufficiently squeezed, it will always remain en- 
tangled independently of the strength of the coupling to 



the environment. The model studied in Ref. 25] is the 
same as our model B but without the interaction be- 
tween the two oscillators. The conclusion in Ref. (2f| 
was derived however using the RWA-Markovian master 
equation. Here we investigate whether the condition pre- 
sented in Ref. [25| is still valid or needs some modification 
in the non-Markovian case. We find that in the case of 
non-Markovian dynamics the condition depends also on 
the interaction strength between the system and environ- 
ment. 

This paper is organized as follows. In Sec. [IT] the 
Hamiltonian and non-Markovian master equations for 
our two models are presented. In Sec. IIIII we introduce 
the concept of renormalization and counter term. In Sec. 
IIVI we introduce the logarithmic negativity to quantify 
entanglement for the two oscillators in our models. This 
entanglement monotone of logarithmic negativity is con- 
veniently computable for general Gaussian states, and 
could provide a proper quantification of entanglement 
particularly for two-mode Gaussian states. We thus dis- 
cuss the covariance matrix and two-mode squeezed vac- 
uum state, a subclass of Gaussian states which we will 
use as initial states. In Sec. El the covariance matrix 
evolution equations obtained from the master equations 
or the Fokker-Planck equations of Wigner function are 
described. In Sec. I VII we present and discuss our re- 
sults of entanglement dynamics based on the logarithmic 
negativity calculated through the evolution equations of 
the covariance matrix. Finally, we investigate whether 
the entanglement survival (or separability time) condi- 
tion under RWA-Markovian approximation in Ref. [2o| 
is still valid in the non-Markovian case. A conclusion is 
given in Sec. I VIII 



II. HAMILTONIAN AND NON-MARKOVIAN 
MASTER EQUATIONS OF TWO MODELS 

In this section we introduce two different but related 
models considered in this paper (see Fig. [T]) and de- 
rive their corresponding quantum master equations up to 
the second order with respect to the system-environment 
coupling constant. The two models both consist of two 
interacting harmonic oscillators. In model A, each of the 
two oscillators is coupled to its own independent ther- 
mal reservoir, while the two oscillators are coupled to a 
common thermal reservoir in model B. 

The Hamiltonian of the system of interest for the two 
models can be written as 

H s = H sl +H s2 + V 12 , (1) 

where H s i and H S 2 are the Hamiltonian of the two sub- 
systems respectively and Viz is the interaction between 
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them. They can be written as 



H S 2 — 
Vv 



Pi 



1 



2M X 2 
Py 1 



MMlx 2 , 



2M y 
Xxy, 



(2) 

(3) 
(4) 



where M x and M y are the masses and £l x and Q y are the 
frequencies of the two subsystems (oscillators), respec- 
tively, and A is the coupling constant between the two 
subsystems. To consider the case of two non-interacting 
oscillators in, we can simply set A to zero. 

We assume that the environments could be described 
as ensembles of harmonic oscillators and interact bilin- 
early through their position operators with the system. 
The Hamiltonian of the two independent environments 
for Model A is thus 

H e = H e \ + H e2 , 

= E(rVr + ^i 1) -i 1)2 ^ )2 ) 



E( Z ^ + ^ 2) ^ ) V 2)2 ), (5) 



and the interaction between the two subsystems and 
reservoirs for Model A is 



V = V1 + V2 



E^^ + EW 



(6) 



where A« and xff are the coupling strengths to their 
own individual reservoir, respectively. 

On the other hand, the Hamiltonian of the environ- 
ment for the two subsystems coupled to a common reser- 
voir (Model B) is 



^=E(|t + ^ 2 a, 



(7) 



and the coupling between two subsystems and reservoirs 
for Model B are 



V = V x +V 2 

= E^^+E^ 2 ^- 



(8) 



Using the perturbative expansion to the second order 
in the system-environment coupling strength, we obtain 
the equation of motion for the reduced density matrix 
p(t) of the system of interest as [35J 



Pit) = -rz[H a ,p{t)\ 
in 



e » 



iH s t 



x(- / dhTr^V^^V^^m^Pe^e^lV) 



where p and V are the density matrix of the system 
and the interaction between the system and environments 
in the interaction picture respectively, and Tr £ indicates 
tracing over environment degrees of freedom with respect 
to the thermal environment density matrix p e . In obtain- 
ing Eq. ((9]), we have also used the fact that the system- 
environment interaction is bilinear in their respective po- 
sitions (or displacements) so that the first-order term 
Tr e (V(t)p e ) = 0. We will derive the non-Markovian 
master equations for our models using Eq. ([9]) without 
making any further approximation. For simplicity, we as- 
sume the masses and resonance frequencies, and coupling 
strengths to the environments are the same for the two 



oscillators. That is M x = M y = M, tt x 



Q y = fi and 



A« = A^ 2 "* = A„. We present the derived non-Markovian 
master equations below. 

Model A: The master equation for two coupled oscil- 
lators, each coupled to its own reservoir, of model A can 
be obtained as 



P(t) 



±-{H s + - A Mtt\{t){x - y) 2 + ~MQl(t){x + y) 2 , p] 
in 4 4 



2h 

i 

1 



Ji(t)[x - y,{p x - p y , p}] 
<y 2 (t)[x + y,{p x +p y ,p}} 



--Di(t)[x-y, [x-y,p]] 
+ ^D 2 (t)[x + y, [x + y,p]} 



2h 
1 

~2h 



fi(t)[x - y, \p x -Py,p]] 
f 2 (t)[x + y, [p x +p y ,p]}, 



(10) 



Here the time-dependent coefficients (lf(t) is called the 
frequency shift, 7i(i) is the dissipation coefficient, and 
Di(t) and fi(t) represent the diffusion coefficients. They 
can be written as 

fi?(t) = ~ J q dt' co S (O i t')'7(* / ) ! (11) 



7<(*) 



1 



Am 



dt'siniriit'W), (12) 



i JO 
t 



Di{t) = - dt' cos(ttit')v(t'), 



(13) 



fi(t) 



1 



MVLi J 



dt'sm(nit')v(t'), (14) 



where i = 1,2, and the frequencies fii and due to the 
interaction A between the two oscillators are 



n 2 



A/M, 



X/M. 



(15) 
(16) 



The two kernels rj(t') and v(t') appearing in Eqs. (jTTJ) 
([14")) are, respectively, the so-called dissipation and noise 
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kernels and are defined as 

1® = ^J2 X n{[Qn(t),q n (0)}) 



dujJ(uj) sin(o>i), 



(17) 



= ^£ A n({<Zn(*),?n(0)}> 



dujJ{uj) cos(wi)(l + 2N(u)), (18) 



where 



N{u>) 



1 



(19) 



is the mean occupation number of the environmental os- 
cillators and 



a^, 



2m n uj n 



S(uj - LU n ) 



(20) 



is the spectral density of the environments. Note that the 
environment position operator q n (t) in Eqs. (|17p and ||18|) 
should be qn\t) for different environments i. We could 
in principle deal with this situation, but for simplicity we 
assume that they have the same corresponding correla- 
tions even though the two environments are independent 
of each other. It is also worth noting that the frequency 
shift (fTTj) and dissipation coefficient (fl~2)) depend only on 
the dissipation kernel (TIT)) while the diffusion coefficients, 
(1131) and (114|) , in the approximation of the perturbative 
expansion depend only on the noise kernel (JTHJ) , thus tem- 
perature dependent. We can see from Eq. (flQ]) that the 
term proportion to ji (t) is responsible for relaxation and 
the term proportional to Di(t) is the main cause for de- 
coherence. The spectral density specifies the structure 
and properties of the environment and thus determines 
the environmental influence on the dynamics of the sys- 
tem of interest. In fact, the time evolution behavior of 
the coefficients of the quantum master equation is rather 
different for environments with different spectral content. 

In principle, we could deal with any given form of the 
spectral density. But as a particular example, we use the 
following form of spectral density to specify the environ- 
ments 136, 37] 



J (to) = — M70CLM 
7T V 



A 



n-l 



2 /A 2 



(21) 



where A is the cutoff frequency, 70 is a constant charac- 
terizing the strength of the interaction with the environ- 
ment, and M is the system mass. The environment is 
said to be ohmic if in the physical range of frequencies 
(u> < A) the spectral density is proportional to uj. And it 
is said to be supra-ohmic if J(u>) is proportional to u) n , 
n > 1, or sub-ohmic if n < 1. For simplicity, in the fol- 
lowing we focus on ohmic baths, i.e., the case of n = 1 in 
Eq. QT). 



If we let A = 0, then this model reduces to two non- 
interacting oscillators, each coupled to its own reservoir. 
The master equation (flQ|) in this instance also reduces to 
the case of just putting two sets of the quantum Brownian 
motion master equations [35| together: 



Pit) 



hH s + lMn 2 (t)(x 2 +y 2 ),p] 
in A 

1 

—p(t)(fa{Px,P}] + [V,{Py,p}]) 
-D{t)([x, [x,p]} + [y, [y,p]}) 

~f(t)([x,\p x ,p]] + [y,lpy,p}]), (22) 



where the time-dependent coefficients f2 2 (£), j{t), D(t), 
and f(t) are defined correspondingly as those in Eqs. (fTT|) - 
HHjl with A = 0. 

Model B: The master equation for two interacting 
oscillators, coupled to a common reservoir, of model B 
can be derived and written as 



U Hs + )-Mhl{t){x + y) 2 ,p] 
in 1 
1 

- jl2 (t) [x + y, {p x + p y , p}} 
-D 2 {t)[x + y, [x + y,p]} 

- T h (t) [X + y, [Px + Py, p}} ■ 



(23) 



The time-dependent coefficients are the same as those in 
Eqs.([Tl~ |) - ([T4]) with i = 2 in model A. Note that com- 
pared with Eq. (HU]), the mode (x — y) is absent in 
Eq. (fl?3"|) . This is a consequence of both the assumptions 
of M x = My = M, Q x = n y = fl and A^ 1} = \i? = A„ 
which we make to simplify the calculation as well as the 
nature of model B which is coupled to a common bath 
with thus the variable (x + y). This can also be inferred 
from an effective factor of two difference in the corre- 
sponding coefficients of terms containing the (x + y) vari- 
able between Eq. (fTOf and Eq. ([23]) . In the derivation of 
the master equation in model B, some sort of additions 
make the coefficients of (x + y) mode twice larger and 
some sort of cancellations make the (x ~ y) mode ab- 
sent in Eq. P3")) . We note that despite being derived 
perturbatively, the master equations for the two mod- 
els seem to be very similar to their exact counterparts 
which are also time-convolutionless with time-dependent 
coefficients 37. M. 



III. RENORMALIZATION AND 
TIME-DEPENDENT COEFFICIENTS 

We note that due to the interaction with the environ- 
ment, the frequency shift term, fl 2 (t) of Eq. (fTTj) . di- 
verges as the cutoff frequency A — > 00 and thus is not 
physical. Thus a regularization procedure is needed for 
the frequency renormalization. There are two different 
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views on this renormalization 39, 40]. One view is start- 
ing from the original Hamiltonian, the frequency can be 
made finite, by a renormalization of frequency, from its 
bare to its physical value. Thus by combining terms which 
involve the frequency in the master equation, the physi- 
cal frequency which is the quantity that can be measured 
in the laboratory is defined as 



n 2 



&(t), 



(24) 



where is the bare frequency in the Hamiltonian. In 
this case, the physical frequency is taken to be finite and 
the bare frequency is taken to be infinite as A — > oo in 
order to cancel the divergent contribution from fl 2 (t). 
Thus the bare frequency has no direct physical signifi- 
cance. Although it may not be exactly the same, this 
view of renormalization has an analogy in solid state 
physics where, for example, electrons are attributed an 
effective mass to take into account their interaction with 
the lattice or/and other electrons. 

An alternative view of renormalization is to regard 
the frequency £1 in the original Hamiltonian as a fi- 
nite renormalized frequency tt r . The fact that this 
Hamiltonian does not give finite frequency then requires 
that extra terms be added to the Hamiltonian to can- 
cel the divergence. These terms are called counter-terms 
[39l l40l l4ll |42| . This view of renormalization is in fact 
more commonly adopted in high energy physics or quan- 
tum field theory. In the context of a quantum Brownian 
motion model, the total Hamiltonian in this case can be 
written as 



H 



1 



(: 



Pn 



^— ' \2m 



1 \ 1 A 2 
+ -m n ujlql + \ n q n x\ + ^ o ( 25 ) 



2 m n tol ' 



1 



Ain 2 r x 2 
2 

n 

2m r- 



X) 



m n LLi 2 



(26) 



The last term in Ec 
counter-term [4lL E 

1 



(125)) can be viewed as a frequency 
with a frequency defined as 



A 



M ^ 



m n uji 



du> 



(27) 



We can find that at large times the frequency shift (l 2 (t) 
is negative and equals to — Q 2 . The added counter-term is 
to cancel the frequency shift at long times and ensure that 
the system can not lower its potential energy below the 
original (renormalized) value. The physical frequency in 
this case is VI 2 = VL 2 +VL 2 +VL 2 (t) and equals to VL 2 at long 
times. We will take this view of renormalization, so the 
original frequency Q in the time-dependent coefficients 
of the master equations, Eqs. (I10|) and (|2"3")) . should be 
replaced by Q — > £l r . 



IV. LOGARITHMIC NEGATIVITY AND 
TWO-MODE GAUSSIAN STATES 

The purpose of this paper is to focus on the entan- 
glement dynamics for the reduced density matrix of the 
two oscillators in our models. The derived master equa- 
tions are generally partial differential equations with 
time-dependent coefficients. Consequently, computing 
the time evolution solution for the density matrix op- 
erator explicitly and then using it to calculate directly 
the dynamics of entanglement are still considered diffi- 
cult. But the problem becomes much more tractable if 
we restrict the states to be Gaussian states. Since the 
couplings in our models are all bilinear in their respective 
positions (displacements) and the effects of the environ- 
ments in the master equations have operator structure no 
more than quadratic in the momenta or/and positions of 
the two oscillators, an initial Gaussian state would re- 
main Gaussian in its subsequent time evolution. So, for 
simplicity, we will consider in the following Gaussian ini- 
tial states for the two oscillators. Any Gaussian state can 
be completely characterized by its corresponding covari- 
ance matrix. We will see below that the time evolution 
of the covariance matrix is easier to calculate than that 
of the density matrix. 

A set of Gaussian states is the set of states with Gaus- 
sian characteristic functions and quasi-probability distri- 
butions of the Wigner function [43[ . The Wigner quasi- 
probability distribution function is defined in terms of 
density matrix p(t) as: 

So a zero-mean Gaussian state is described, for example, 
by the Wigner function as 



W(K) 



- ~ ) A = expt-lxV^X 7 ), (29) 
47r n VdetV 2 



where V is the covariance matrix and X represents the 
vector {x\,pi, x%, £>2> •■ • iX n ,p n ). The zero mean denotes 
(Xi) = 0, and this can be changed at will using local 
unitary displacement operators. So we can set (JQ) = 
without loss of generality. 

The matrix elements of the covariance matrix V are 
defined as 

Vij = ({AX^AXj}) =Tr({AA > l ,AA > J } ( 5) 

= J d 4 X AXiAXjW(X), (30) 

where {AX^AXj} = (AXiAXj + AX J AX l )/2, AX, = 
Xi — (Xi) and AXi = Xi — (Xi). The average of the 
operator Xi, (Xi), means Tr(X iy o), and (Xi) denotes an 
average of a variable Xi with respect to the Wigner func- 
tion distribution W(X), so (Xi) equals to (Xi). With 
this definition, we could transfer the problem of solving 
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a time-dependent partial differential equation of the den- 
sity matrix into a problem of solving first-order in time, 
coupled linear ordinary differential equations of the co- 
variance matrix elements. This could be done by first 
transferring the master equations to the Fokker-Planck 
equation for the Wigner function, and then finding the 
coupled differential evolution equation for the covariance 
matrix elements using Eq. I|30p and the Fokker-Planck 
equation (see Sec. |V|) . 

We will use the logarithmic negativity to quantify the 
degrees of entanglement of the infinite-dimensional bipar- 
tite system states of the two oscillators. The logarithmic 
negativity of a bi par tite system consisting of two subsys- 
tems A and B is |44| 



where 



E N (p) = log 2 || p 



(31) 



where p Ts means partial transpose of a (mixed) state 
density matrix operator p with respect to subsystem B. 
That is to say, (i A , ]b\p Tb |&a ,1b) = Jb) for 

any arbitrary orthonormal product basis which is be- 
longed to the tensor product of Hilbert space of com- 
binative system A and B. The operation || . ||i denotes 
the trace norm and the trace norm of any Hermitian op- 
erator H is defined as || H ||i= Tr|if| = TtVhTH. 

Despite not being convex, the logarithmic negativity 
is a full entanglement monotone under local operations 
and classical communication [45| and constitutes to an 
upper bound to the distillable entanglement [44j| . For 
the particular case of two-mode Gaussian states, the log- 
arithmic negativity could actually provide an appropri- 
ate quantification of quantum entanglement. Vidal and 
Werner [44[ demonstrated that logarithmic negativity is 
computable for general Gaussian states. For two-mode 
Gaussian states, it can be furthermore shown that the 
logarithmic negativity can be represented as [43[ 



E N (p) = max(0, - log 2 2V S 



(32) 



where V s is the smallest sympletic eigenvalue of the 
partially transposed covariance matrix of the two-mode 
Gaussian states. Equation (f3"2"| is a simple decreasing 
function of V s which quantifies the degree of violation of 
the necessary and sufficient separability criterion of the 
positivity of partial transpose [431 , H3| ■ For > 1/2 the 
state is separable, otherwise it is entangled. So the small- 
est partially transposed sympletic eigenvalue V s alone 
completely qualifies and quantifies the quantum entan- 
glement of a two-mode Gaussian state [43j |. That is, the 
smaller the value of V s , the more entangled the corre- 
sponding two-mode Gaussian state. As a result, the log- 
arithmic negativity may be regarded as a suitable entan- 
glement quantification indicator for two-mode Gaussian 
states. 

The partially transposed sympletic eigenvalues Vi are 
the symplectic eigenvalues of V Tfi , and V Ts can be writ- 
ten down in a compact form [l8j 



P 



1 
1 



f 

-f 



(34) 



and A®B means that the block-diagonal matrix with the 
matrices A and B as diagonal entries. The symplectic 
eigenvalues are the positive square roots of the standard 
eigenvalues of the —a\ TB uY Tb or the absolute value of 
the eigenvalues of ioY TB . Here a is called the symplectic 
matrix from the commutation relations = ihuij 

which is given by 



J 
J 



and J 



f 

-f 



(35) 



The logarithmic negativity in the form of Eq. (f3"2"| is 
much more easier to compute than that defined in Eq. 
(EED- 

One subclass of two-mode Gaussian states is the so- 
called two-mode squeezed vacuum states. The posi- 
tion and momentum wave functions for the two-mode 
squeezed vacuum state with a squeezing parameter r are 



.1] 



i>(x,y) = ;/~exp[-e 2r (x + yf/2 



1p(j>x,Py) 



e 2 ^(x-y) 2 /2}, 
2 

exp[~e- 2r ( Px -pyf/2 
-e 2r (p x+Pv ) 2 /2} 



(36) 



(37) 



They approach C6(x — y) and C5(p x +p y ) , respectively, 
in the limit of infinitely squeezing r — > oo, where C is 
some constant. The corresponding Wigner function of 
the two-mode squeezed vacuum state is then [l[ 



W(X) = —eM^- 2r [(x + y) 2 + (p x - Py ) 2 

7T 

- e ^( X -yf + (p x+ p y ) 2 ]}. 



(38) 



(33) 



In the limit of infinitely squeezing r — > oo, this Wigner 
function approaches CS(x — y)5(p x + p y ), corresponding 
to the original (perfectly correlated and maximally en- 
tangled) EPR state. While at r = 0, the two-mode state 
corresponds to a separable (disentangled) state. The two- 
mode squeezed vacuum states are routinely generated in 
quantum optics laboratories and have been used in most 
implementations of continuous variable quantum infor- 
mation protocols P, 0, 5f|. It has also been proposed 
recently that a two-mode squeezed state could be gener- 
ated for two nanomechanical oscillators that act as the 
two opposite sections, suspended from the substrate, of 
a dc-SQUID (superconducting quantum interference de- 
vice) loop [la ]. 

The two-mode squeezed vacuum states, from Eqs. (|2T))) 
and (|3"5|) , can be completely characterized by the follow- 
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ing covariance matrix [47|, l48| : 



the renormalized frequency Q — > fl r 



(a -c 

a c 

-c a 

V c a. 



(39) 



where a and c are 



cosh(2r) /2, 



sinh(2r)/2. 



(40) 



For simplicity, we will use the two-mode squeezed vacuum 
states as the initial states of the two quantum oscillators 
in our models throughout the paper. 



V. EVOLUTION EQUATIONS OF THE 
COVARIANCE MATRIX ELEMENTS 



V X1 = 

Vl2 = 

Vl3 = 

V u = 

v 22 = 

v 23 = 



V: 



21 



^33 
V34 



2V 12 , 

-(Q 2 + Q 2 C + n%{t))V xl - 2 l2 {t)V l2 

-(n 2 c + n|(t) + x/M)Via - 2 l2 {t)v u + V22 - hf 2 (t), 

Vu + V 23 , 

-(Q 2 C + n\{t) + X/M)Vu - 2 l2 {t)V l2 

-(n 2 r + n 2 c + fi 2 2 (t))v 13 - 2 l2 (t)v u + v 2i - hf 2 (t), 



-ml 
-ml 
-(n?-f 
-(VI- 



Vm = 



+ n 2 c + n 2 2 (t))v 12 - 4 72 (t)v 22 

F n 2 (t) + X/M)V 23 - ij 2 (t)V 24 + 2H 2 D 2 (t), 

n 2 c + n 2 (t))v 13 - 2 l2 {t)v 23 + v 2i 

& 2 + X)V 33 -2 l2 V 3i -hf 2 (t), 

-{n 2 c + nj(t) + x/M)v u - (n 2 + n 2 c + n§(t))yu 
-2 l2 {t)v 22 - (n 2 r + n 2 c + Ci 2 {t))v 2 3 - 4 7 2(t)F 24 
-(nl + n 2 {t) + x/m)v 34 - 2 72 (i)V44 + 2h 2 D 2 {t), 

2V 34l 

-(n 2 c + n 2 (t))v 13 -2 l2 (t)v 23 

-(Q 2 r + n 2 c + n 2 2 (t))v 33 - 2 l2 (t)v 3i + f 44 - hf 2 (t), 

-2(n 2 c + n 2 2 (t) + x/M)v u - 4 72 (t)V24 

-2(n? + n 2 c + n!(t))Vs4 - 4:j 2 (t)v u + 2h 2 D 2 (t). 



Similar calculations can be performed for model A. The 
solutions of the coupled first-order ordinary differential 
equations in time are much easier to calculate than the 
partial differential equations of the Fokker-Planck equa- 
tions or the quantum master equations. Solving for the 
time evolution of the covariance matrix elements, we can 
then obtain the entanglement dynamics through the com- 
putation of the logarithmic negativity using Eq. (f32|) . 



By the definition of Wigner function in Eq. the 
corresponding function for p(t) is 



^^*)=(i)7I d ^-f | ^ )l9+ f )exp( ^ 1) - 

(41) 

Using Eqs. f28]) and (f4Tj) , and after straightforward but 
somehow tedious calculations, we can obtain the Fokker- 
Planck equations of the Wigner function corresponding 
to the master equations (fit))) and (j2"3")) . From the Fokker- 
Planck equation of the Wigner function and Eq. ([30]) , we 
can obtain coupled first-order ordinary differential equa- 
tions with time-dependent coefficients for all elements of 
the covariance matrix. Due to the symmetrical property 
of the covariance matrix, i.e. V\ 2 — V 2 i, V24 = V42 etc., 
we only need ten essential components of the covariance 
matrix in the bipartite system here instead of sixteen 
components. For example, we obtain for model B with 



VI. RESULTS AND DISCUSSIONS 

We first report our numerical results of the non- 
Markovian entanglement (logarithmic negativity) dy- 
namics for our two models in two cases: (1) 70 = (iso- 
lated), (2) 70 = 6 x 10 _2 r2 r , where 70 is a constant in 
the spectral density (f2"Tj) and is related to the coupling 
strength to the environments and f2 r is the renormalized 
frequency of the subsystems. In all the plots presented 
below, the parameters used are as follows. The environ- 
ment temperature is at kgT = 10HQ r , the cutoff fre- 
quency is A = 2000J7 r and the interaction between two 
subsystems A is in units of Mfl 2 ,. Finally, we investi- 
gate whether the entanglement survival condition under 
RWA-Markovian approximation in Ref. [25( is still valid 
in the non-Markovian case. 



A. Isolated system (70 = 0) 

Similar calculations can be performed for model A. 
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FIG. 2: (Color online). Time evolution of the logarithmic 
negativity of the two subsystems isolated from external envi- 
ronments (70 = 0) for the case of an initial two-mode squeezed 
vacuum state with a squeezing parameter r = 2. The solid 
line stands for A = 0, dashed for A = 0.2, dotted for A = 0.8, 
dash-dotted for A = —0.2 and dash-dot-dotted for A = —0.8. 
In all the plots presented below, the parameters used are as 
follows. The environment temperature is at fcsT = 10M7 r , 
the cutoff frequency is A = 2000f2 r and the interaction be- 
tween two subsystems A is in units of MfiJL 
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FIG. 3: (Color online). Time evolution of the logarithmic 
negativity for an initial two-mode squeezed vacuum state with 
a squeezing parameter r — 0.1. Other conditions and plot 
caption are the same as in Fig. [2] 



We plot the dynamics of the logarithmic negativity of 
the two models when they are isolated from the external 
environments in Figs. [2HU Similar to that of an inter- 
acting discrete system of two qubits, the dynamics of 
entanglement of the two oscillators depends strongly on 
the initial states and on the interacting strength between 



FIG. 4: (Color online). Time evolution of the logarithmic neg- 
ativity for an initial two-mode squeezed vacuum state with a 
squeezing parameter r = 0. Other conditions and plot caption 
are the same as in Fig. 

them. When there is no interaction between the two 
subsystems (A = 0) isolated from the environments, the 
time evolution of the logarithmic negativity maintains 
constant in Figs. HHH as it should. But the dynamics 
of the logarithmic negativity varies quasi-periodically for 
two interacting subsystems isolated from the external en- 
vironments and the smaller the value of the interaction 
strength |A| between the two oscillators is, the longer the 
quasi-period of the logarithmic negativity is. This can 
be seen from the plots of the entanglement dynamics in 
a longer time scale. 

From Fig. 01 we see that the entanglement of the two 
subsystems can be generated from an initially separa- 
ble state (r = 0) through their mutual interaction, and 
the larger the interaction strength, the larger the gener- 
ation of the entanglement. We can also see from Fig. 0] 
that the entanglement dynamics for an initially separa- 
ble state (r = 0) seems to be symmetric with respect to 
the change of Xxy <-> —Xxy. While this is not the case 
for r ( see Figs. [2] and [3]). This may be due to the 
fact that for r — the initial wave function, Eq. (|36p . 
or the Wigner function, Eq. (|38[) . is symmetric under the 
change of x <-> —x, y <-> — y and xy <-> — xy. While 
for r ^ 0, this symmetry is broken and the initial wave 
function, Eq. (f3"rj)) . or the Wigner function, Eq. ((3"5)) . pos- 
sesses the preferred entanglement in the relative position 
variable (x — y) as compared to the variable (x + y). So 
an attractive interaction (A < 0) seems to enhance this 
entanglement in (x — y). This can be seen from Figs. [2HH 
that for a fixed value of the squeezing parameter r, if the 
interaction strength is attractive (A < 0) then the entan- 
glement grows initially with time. On the other hand, the 
entanglement decreases with time initially if the interac- 
tion strength is positive and smaller than a critical value, 
i.e. < A < A c . For example, in Fig. [3] the initial entan- 
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FIG. 5: The positive interaction strength versus initial squeez- 
ing parameter phase diagram. The regime above (below) the 
critical line curve corresponds to the situation that the entan- 
glement increases (decreases) with time initially. 



glement grows with time for A = 0.8 while it decreases 
with time for A = 0.2. Similarly, we may say that for a 
fixed positive value of A > 0, there exists a critical ini- 
tial squeezing parameter above which the entanglement 
decreases with time initially. Figure [5] shows the criti- 
cal value line that separates these two situations in the 
positive interaction strength A versus initial squeezing 
parameter r phase diagram. 

The general trend is that when the entanglement grows 
with time initially, the entanglement is enhanced to reach 
maximum values at later times; while if the entanglement 
decreases with time initially, the initial value of the en- 
tanglement is usually the maximum value. Thus, for two 
oscillators initially in a two-mode squeezed vacuum state 
with a squeezing parameter r, an attractive interaction 
(A < 0) between the two oscillators is able to enhance 
their entanglement at later times. An repulsive inter- 
action (A > 0) can enhances the entanglement at later 
times if A > A c , but the entanglement is no longer in- 
creased if < A < A c . For the same value of interaction 
strength between the two oscillators, the attractive in- 
teraction seems always to be better than the repulsive 
interaction as far as the maximum value of entanglement 
that can be reached at a later time is concerned. 



B. Coupled to environments (70 = 6 x 10 2 Cl r ) 

In Figs. [6] and [3 we plot the dynamics of the loga- 
rithmic negativity of the two subsystems coupled to the 
environments for our two models with different initial 
states of r — 2 and r = 0, respectively. Compared with 
the corresponding 70 = cases in Figs. [2] the oscillatory 
phenomena due to the influence of the environments dis- 
appear except in Fig. Efb). It can also be seen from Fig. 
[SJa) that the entanglement vanishes in finite times (sud- 
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FIG. 6: (Color online). Time evolutions of the logarithmic 
negativity of the two subsystems coupled more strongly to 
the environments (70 = 6x 10 _2 f2 r ) for an initial two-mode 
squeezed vacuum state with a squeezing parameter r = 2. 
The subplot (a) is for Model A and (b) is for Model B. The 
solid line is for A = 0, dashed for A = 0.2, dotted for A = 0.8, 
dash-dotted for A = —0.2 and dash-dot-dotted for A = —0.8. 



den death) [49( . This is in contrast to the loss of quantum 
coherence that is usually gradual [3(1 HH H§| . 

The logarithmic negativity, shown in Fig. [BJa) decays 
very fast for model A as compared to Model B in Fig. 
[B^b). In other words, the entanglement can sustain much 
longer when two subsystems are coupled to a common 
bath than to individually independent baths. This con- 
clusion is consistent with the result found in other con- 
tinuous variable models [2^, [28| or discrete qubit models 
[HTI . [52I ]. In our models, this could be understood by 
noting that the Hamiltonian of the total system can be 
written in terms of new dynamical variables, the sum and 
difference of the two oscillator's positions and momenta 
(x + y, p x + p y , x — y, p x — p y ). For model B coupled 
to a common environment, only the mode of the sum of 
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FIG. 7: (Color online). Time evolutions of the logarithmic 
negativity for an initial squeezing parameter r = 0. Other 
conditions and plot caption are the same as in Fig. [6] 



the two positions interacts with the environment and the 
mode of the difference of the two positions undergoes a 
free evolution. As a result, only the modes of the sum of 
the positions and momenta are affected by the environ- 
ment [see Eqs. (|2"3")) ]. But for Model A, these modes all 
interact with the environments [see Eqs. (|10p] and thus 
are all influenced by the environments. 

In Fig. [7fa), we find that the logarithmic negativity 
for Model A is barely generated with an initially sepa- 
rable state (r = case). On the other hand, we find 
in Fig. [7fb) that even with no interaction between the 
two subsystems, but due to the fact that they coupled 
to a common bath (Model B), the entanglement can be 
generated from an initially separable state (r = case) 
pH HI, [51], [52| . But the generated entanglement lasts 
only for a short time and then disappears. In most situ- 
ations the entanglement is created for a very short time 
after the interaction with the environment is turned on. 
The entanglement may persist for long times or disap- 
pear shortly, depending on system-environment coupling 
and the properties of the environment [III, [28|, [5l], ■ 



C. Condition for Entanglement Survival 

A condition derived in Ref. [25| stated that if the 
two-mode squeezed state of the electromagnetic field em- 
bedded in a thermal environment is initially sufficiently 
squeezed, it will always remain entangled independently 
of the strength of the interaction to the environment. 
Each of the two electromagnetic modes has the Hilbert 
space equivalent to the Hilbert space of the Harmonic os- 
cillator. As a result, the model studied in Ref. [25| is the 
same as our Model B during A = 0. However, the conclu- 
sion in Ref. (2f| was reached using the RWA-Markovian 
master equation. Here we investigate whether the con- 
dition presented in Ref. [25[ is still valid or needs some 
modification in the non-Markovian case. 

In Ref. [25|, the Simon criterion [47| for continuous 
variables system was used to verify whether the quantum 
state of the system is entangled or separable. It was 
found that if the initial state is sufficiently squeezed, i.e., 
the squeezed parameter of the initial two-mode squeezed 
vacuum state satisfies 12511. 



|r| ^ -ln(2A + l), 



where 



N 



1 



(42) 



(43) 



is a mean thermal photon number, it will remain entan- 
gled forever in spite of the interaction between the system 
and the external environment. Otherwise, the state will 
disentangle (become separable) after time |25[ 



1 , /2N + l-e 
— In — = 

2 7 V 2N 



-2\r\ 



1 



2\r\ 



(44) 



where 7 w 21im t _ >00 7(f) w 270 From Eqs. (|42|) and ([43]) 
at a temperature of ksT = 10M7 r , the corresponding 
critical squeezed parameter is \r c \ = iln(2iV + 1) = 
1.498. We choose the squeezing parameters to be at 
and slightly smaller than this critical squeezing value, 
and vary the system-environment interaction strengths 
to check whether the condition for the inequality (|42]) is 
still valid. 

Form Fig. [51(a) and its inset the logarithmic negativity 
does not vanish and appear cyclically at long time for 
r = r c = 1.498 regardless of their system-environment 
interaction strength. So the statement about the inequal- 
ity (|4"2"| in Ref. [25[ seems valid for both non-Markovian 
and RWA-Markovian cases. On the other hand, if the 
squeezed parameter r = 1.4 < r c , i.e., smaller than 
the critical squeezed parameter, Eq. (|4"4")l predicts that 
the two-mode state will disentangled (or become sepa- 
rable) after time t = 7.168/f2 r for 70 = 0.06f2 r and 
t = 0.43/O r for 70 = iri r , respectively. This is indeed the 
case for RWA-Markovian approximation results shown in 
Fig. [5Jb) for our Model B with A = 0. However, this is 
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FIG. 8: (Color online). The time evolutions of the logarith- 
mic negativity for different initial squeezed parameters, (a) 
r=1.489 and (b) r=1.4. Two different values of 70 (70 = 
0.06fi r and 70 = l£l r ) are used in each plot, where 70 is 
related to the system-environment coupling strength. The in- 
sets illustrate the same plots but with much smaller values of 
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system and the environment. 

VII. CONCLUSION 

We have investigated the non-Markovian entanglement 
dynamics of two oscillator subsystems which are coupled 
to a common environment or are coupled respectively 
to their own independent environments. We have pre- 
sented and discussed the influence of the environments 
on the entanglement dynamics by varying initial states 
(different squeezing parameters), oscillator-oscillator in- 
teractions and oscillator-environment interactions. We 
have found that the dynamics of entanglement oscillates 
for two interacting subsystems isolated from the external 
environments. The attractive interaction seems always 
to be better than the repulsive interaction as far as the 
maximum value of entanglement that can be reached at 
a later time is concerned. When the coupling between 
the environments and the two subsystems is turned on 
and increased progressively, these periodic behaviors die 
down gradually and disappear eventually. When the in- 
teraction strength to the environments is increased fur- 
ther, the entanglement vanishes at finite times (sudden 
death). This is in contrast to the loss of quantum coher- 
ence that is usually gradual. It is also been found that 
the entanglement can sustain much longer when the two 
subsystems are coupled to a common bath than to indi- 
vidually independent baths. This conclusion is consistent 
with the result found in other models [25| . In summary, 
the dynamics of the quantum entanglement is sensitive 
to the initial states, the oscillator-oscillator interaction, 
the oscillator-environment interaction and the coupling 
to a common bath or to different, independent baths. 

Finally, we have checked the condition for entangle- 
ment to exist at long times for two non-interacting sub- 
systems coupled to a common bath (model B with A = 0). 
In contrast to the condition, which depends only on the 
mean thermal phonon number, obtained in Ref. [2 a ] using 
RWA-Markovian master equation, our non-Markovian 
analysis indicates that the condition also depends on the 
system-environment interaction. 



not true for the non-Markovian case. We find that en- 
tanglement disappear except for the non-Markovian case 
with a larger coupling strength 70, in which the entangle- 
ment dies out firstly and then be generated cyclically by 
the interaction to the common bath [see Fig. Hfb) and 
its inset]. In other words, the non-Markovian dynamics 
predicts that the entanglement would persist for a longer 
time. This is consistent with the result in Ref. [Hj]. So 
in the case of non-Markovian dynamics, the inequality 
(p4"2")) and Eq. are no longer true, and the condition 
not only depends on the mean thermal photon number 
but also depends on the interaction strength between the 
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